Intro

The first example is about Geometric people. This is an excuse to revisit models in which the relationship among parameters is not necessarily linear. One of the main benefits is that these models can be constructed following theoretical concepts and thus, the parameters have meaningful interpretations. There are somethings to keep in mind:

Parameters need to be identifiable. If this is not obviously mathematically, it will be when you try to fit the model. The parameters will be correlated to a high degree or the model will not converge. In the HMC context, you will get divergence or other small ESS or something else.
You might have competing models which are very reasonable a priori and will need to fit them and compare them using WAIC or other criteria.
You can assign priors for these type of models based on previous data or experience. You do not have to always scale variables as it is done in the book

Nonlinear models: Simple example using an asymptotic regression model

data(barley, package = "nlraa")

## Visualize
ggplot(data = barley, aes(x = NF, y = yield)) + 
  geom_point() + ylab("Yield (g/m2)") + xlab("Nitrogen Fertilizer (g/m2)")

priors <- prior(normal(400, 100), nlpar = "Asym", coef = "Intercept") + 
  prior(normal(-2, 2), nlpar = "lrc", coef = "Intercept") + 
  prior(normal(100, 50), nlpar = "R0", coef = "Intercept")

bf1 <- bf(yield ~ Asym + (R0 - Asym)*exp(-exp(lrc) * NF), 
           Asym + R0 + lrc ~ 1,
           nl = TRUE)

brm1 <- brm(bf1, data = barley, seed = 123, prior = priors, refresh = 0)

## Compiling Stan program...

## Start sampling

### Output
brm1

##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: yield ~ Asym + (R0 - Asym) * exp(-exp(lrc) * NF) 
##          Asym ~ 1
##          R0 ~ 1
##          lrc ~ 1
##    Data: barley (Number of observations: 76) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Asym_Intercept   395.89     41.24   337.24   497.41 1.00     1614     1017
## R0_Intercept     132.33     15.69   102.00   163.19 1.00     2414     2153
## lrc_Intercept     -1.78      0.36    -2.52    -1.10 1.00     1538      993
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma    72.26      6.23    61.65    86.16 1.00     2727     2030
## 
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

## Visualize results
plot(brm1, "^b_")

## Pairs plot
pairs(brm1, "^b_")

## Can we get the correlation matrix for this plot?
round(cov2cor(vcov(brm1)), 2)

##                Asym_Intercept R0_Intercept lrc_Intercept
## Asym_Intercept           1.00         0.17         -0.90
## R0_Intercept             0.17         1.00         -0.31
## lrc_Intercept           -0.90        -0.31          1.00

## Plot conditinal effects
plot(conditional_effects(brm1), points = TRUE)

## How do I just extract the parameter estimates?
prm.est <- fixef(brm1)
round(prm.est, 1)

##                Estimate Est.Error  Q2.5 Q97.5
## Asym_Intercept    395.9      41.2 337.2 497.4
## R0_Intercept      132.3      15.7 102.0 163.2
## lrc_Intercept      -1.8       0.4  -2.5  -1.1

## There is a Bayes R^2, which I guess it means something??
bayes_R2(brm1)

##     Estimate  Est.Error      Q2.5    Q97.5
## R2 0.6023969 0.04468204 0.5016564 0.672277

## Compute the WAIC and LOO
waic(brm1)

## Warning: 
## 1 (1.3%) p_waic estimates greater than 0.4. We recommend trying loo instead.

## 
## Computed from 4000 by 76 log-likelihood matrix
## 
##           Estimate   SE
## elpd_waic   -434.4  7.1
## p_waic         3.7  0.8
## waic         868.9 14.1
## 
## 1 (1.3%) p_waic estimates greater than 0.4. We recommend trying loo instead.

loo(brm1)

## 
## Computed from 4000 by 76 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo   -434.5  7.1
## p_loo         3.7  0.8
## looic       868.9 14.2
## ------
## Monte Carlo SE of elpd_loo is 0.0.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.

Hidden Minds and Observed Behavior

The so-called inverse problem is one of the most basic problems in scientific inference: How to figure out causes from observations. It is a problem, because many different causes can produce the same evidence.

Boxes example

data(Boxes_model)
cat(Boxes_model)

## 
## data{
##     int N;
##     int y[N];
##     int majority_first[N];
## }
## parameters{
##     simplex[5] p;
## }
## model{
##     vector[5] phi;
##     
##     // prior
##     p ~ dirichlet( rep_vector(4,5) );
##     
##     // probability of data
##     for ( i in 1:N ) {
##         if ( y[i]==2 ) phi[1]=1; else phi[1]=0; // majority
##         if ( y[i]==3 ) phi[2]=1; else phi[2]=0; // minority
##         if ( y[i]==1 ) phi[3]=1; else phi[3]=0; // maverick
##         phi[4]=1.0/3.0;                         // random
##         if ( majority_first[i]==1 )             // follow first
##             if ( y[i]==2 ) phi[5]=1; else phi[5]=0;
##         else
##             if ( y[i]==3 ) phi[5]=1; else phi[5]=0;
##         
##         // compute log( p_s * Pr(y_i|s )
##         for ( j in 1:5 ) phi[j] = log(p[j]) + log(phi[j]);
##         // compute average log-probability of y_i
##         target += log_sum_exp( phi );
##     }
## }

Ordinary Differential Equations

For some background: https://en.wikipedia.org/wiki/Ordinary_differential_equation

Nut cracking

The simple equation which describes nut cracking is

\[ \frac{dM}{dt} = k (M_{max} - M_t) \]

Using rethinking

data("Panda_nuts")

## R code 16.11
dat_list <- list(
    n = as.integer( Panda_nuts$nuts_opened ),
    age = Panda_nuts$age / max(Panda_nuts$age),
    seconds = Panda_nuts$seconds )

m16.4 <- ulam(
    alist(
        n ~ poisson( lambda ),
        lambda <- seconds*phi*(1-exp(-k*age))^theta,
        phi ~ lognormal( log(1) , 0.1 ),
        k ~ lognormal( log(2) , 0.25 ),
        theta ~ lognormal( log(5) , 0.25 )
    ), data=dat_list , chains=4 )

## 
## SAMPLING FOR MODEL '36b13a4a3d6bae56fb4ea916b4826328' NOW (CHAIN 1).
## Chain 1: 
## Chain 1: Gradient evaluation took 5.1e-05 seconds
## Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.51 seconds.
## Chain 1: Adjust your expectations accordingly!
## Chain 1: 
## Chain 1: 
## Chain 1: Iteration:   1 / 1000 [  0%]  (Warmup)
## Chain 1: Iteration: 100 / 1000 [ 10%]  (Warmup)
## Chain 1: Iteration: 200 / 1000 [ 20%]  (Warmup)
## Chain 1: Iteration: 300 / 1000 [ 30%]  (Warmup)
## Chain 1: Iteration: 400 / 1000 [ 40%]  (Warmup)
## Chain 1: Iteration: 500 / 1000 [ 50%]  (Warmup)
## Chain 1: Iteration: 501 / 1000 [ 50%]  (Sampling)
## Chain 1: Iteration: 600 / 1000 [ 60%]  (Sampling)
## Chain 1: Iteration: 700 / 1000 [ 70%]  (Sampling)
## Chain 1: Iteration: 800 / 1000 [ 80%]  (Sampling)
## Chain 1: Iteration: 900 / 1000 [ 90%]  (Sampling)
## Chain 1: Iteration: 1000 / 1000 [100%]  (Sampling)
## Chain 1: 
## Chain 1:  Elapsed Time: 0.158473 seconds (Warm-up)
## Chain 1:                0.137153 seconds (Sampling)
## Chain 1:                0.295626 seconds (Total)
## Chain 1: 
## 
## SAMPLING FOR MODEL '36b13a4a3d6bae56fb4ea916b4826328' NOW (CHAIN 2).
## Chain 2: 
## Chain 2: Gradient evaluation took 1.5e-05 seconds
## Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.15 seconds.
## Chain 2: Adjust your expectations accordingly!
## Chain 2: 
## Chain 2: 
## Chain 2: Iteration:   1 / 1000 [  0%]  (Warmup)
## Chain 2: Iteration: 100 / 1000 [ 10%]  (Warmup)
## Chain 2: Iteration: 200 / 1000 [ 20%]  (Warmup)
## Chain 2: Iteration: 300 / 1000 [ 30%]  (Warmup)
## Chain 2: Iteration: 400 / 1000 [ 40%]  (Warmup)
## Chain 2: Iteration: 500 / 1000 [ 50%]  (Warmup)
## Chain 2: Iteration: 501 / 1000 [ 50%]  (Sampling)
## Chain 2: Iteration: 600 / 1000 [ 60%]  (Sampling)
## Chain 2: Iteration: 700 / 1000 [ 70%]  (Sampling)
## Chain 2: Iteration: 800 / 1000 [ 80%]  (Sampling)
## Chain 2: Iteration: 900 / 1000 [ 90%]  (Sampling)
## Chain 2: Iteration: 1000 / 1000 [100%]  (Sampling)
## Chain 2: 
## Chain 2:  Elapsed Time: 0.164673 seconds (Warm-up)
## Chain 2:                0.126106 seconds (Sampling)
## Chain 2:                0.290779 seconds (Total)
## Chain 2: 
## 
## SAMPLING FOR MODEL '36b13a4a3d6bae56fb4ea916b4826328' NOW (CHAIN 3).
## Chain 3: 
## Chain 3: Gradient evaluation took 1.5e-05 seconds
## Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.15 seconds.
## Chain 3: Adjust your expectations accordingly!
## Chain 3: 
## Chain 3: 
## Chain 3: Iteration:   1 / 1000 [  0%]  (Warmup)
## Chain 3: Iteration: 100 / 1000 [ 10%]  (Warmup)
## Chain 3: Iteration: 200 / 1000 [ 20%]  (Warmup)
## Chain 3: Iteration: 300 / 1000 [ 30%]  (Warmup)
## Chain 3: Iteration: 400 / 1000 [ 40%]  (Warmup)
## Chain 3: Iteration: 500 / 1000 [ 50%]  (Warmup)
## Chain 3: Iteration: 501 / 1000 [ 50%]  (Sampling)
## Chain 3: Iteration: 600 / 1000 [ 60%]  (Sampling)
## Chain 3: Iteration: 700 / 1000 [ 70%]  (Sampling)
## Chain 3: Iteration: 800 / 1000 [ 80%]  (Sampling)
## Chain 3: Iteration: 900 / 1000 [ 90%]  (Sampling)
## Chain 3: Iteration: 1000 / 1000 [100%]  (Sampling)
## Chain 3: 
## Chain 3:  Elapsed Time: 0.14697 seconds (Warm-up)
## Chain 3:                0.150362 seconds (Sampling)
## Chain 3:                0.297332 seconds (Total)
## Chain 3: 
## 
## SAMPLING FOR MODEL '36b13a4a3d6bae56fb4ea916b4826328' NOW (CHAIN 4).
## Chain 4: 
## Chain 4: Gradient evaluation took 1.7e-05 seconds
## Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.17 seconds.
## Chain 4: Adjust your expectations accordingly!
## Chain 4: 
## Chain 4: 
## Chain 4: Iteration:   1 / 1000 [  0%]  (Warmup)
## Chain 4: Iteration: 100 / 1000 [ 10%]  (Warmup)
## Chain 4: Iteration: 200 / 1000 [ 20%]  (Warmup)
## Chain 4: Iteration: 300 / 1000 [ 30%]  (Warmup)
## Chain 4: Iteration: 400 / 1000 [ 40%]  (Warmup)
## Chain 4: Iteration: 500 / 1000 [ 50%]  (Warmup)
## Chain 4: Iteration: 501 / 1000 [ 50%]  (Sampling)
## Chain 4: Iteration: 600 / 1000 [ 60%]  (Sampling)
## Chain 4: Iteration: 700 / 1000 [ 70%]  (Sampling)
## Chain 4: Iteration: 800 / 1000 [ 80%]  (Sampling)
## Chain 4: Iteration: 900 / 1000 [ 90%]  (Sampling)
## Chain 4: Iteration: 1000 / 1000 [100%]  (Sampling)
## Chain 4: 
## Chain 4:  Elapsed Time: 0.142407 seconds (Warm-up)
## Chain 4:                0.118888 seconds (Sampling)
## Chain 4:                0.261295 seconds (Total)
## Chain 4:

## What are the parameter estimates?
precis(m16.4)

##            mean         sd      5.5%      94.5%    n_eff    Rhat4
## phi   0.8662389 0.03931907 0.8043098  0.9282628 789.8238 1.002701
## k     5.9633572 0.56082312 5.0613425  6.8575892 601.9944 1.003945
## theta 9.8261385 2.02504422 6.9313479 13.3317491 620.4840 1.004014

Using brms

Panda_nuts$s_age <- Panda_nuts$age / max(Panda_nuts$age)
## Setting up priors
nut_prs <- prior(lognormal(log(1), 0.1), nlpar = "phi", coef = "Intercept") + 
  prior(lognormal(log(2), 0.25), nlpar = "k", coef = "Intercept") + 
  prior(lognormal(log(5), 0.25), nlpar = "theta", coef = "Intercept")

nut_bf <- bf(nuts_opened ~ seconds * phi * (1 - exp(-k * s_age))^theta,
             phi + k + theta ~ 1, nl = TRUE, family = poisson(link = "identity"))

nuts_brm <- brm(nut_bf,
                data = Panda_nuts,
                seed = 123,
                prior = nut_prs,
                refresh = 0)

## Warning: It appears as if you have specified a lower bounded prior on a parameter that has no natural lower bound.
## If this is really what you want, please specify argument 'lb' of 'set_prior' appropriately.
## Warning occurred for prior 
## b_phi_Intercept ~ lognormal(log(1), 0.1)
## b_k_Intercept ~ lognormal(log(2), 0.25)
## b_theta_Intercept ~ lognormal(log(5), 0.25)

## Compiling Stan program...

## Start sampling

## Looks like now we get the same answer
plot(nuts_brm, "^b_")

## Pairs plot
pairs(nuts_brm, "^b_")

## How does the plot look?
plot(conditional_effects(nuts_brm), ask = FALSE, points = TRUE)

## Predictions
# ndat <- expand.grid(s_age = seq(0, 1, length.out = 20), seconds = c(2, 20, 40, 80, 130))
# prds <- predict(nuts_brm, newdata = ndat)
# prds2 <- merge(ndat, prds)
# prds2$seconds_f <- as.factor(prds2$seconds)
# prds2$age <- prds2$s_age * max(Panda_nuts$age)
# ## Plot
# ggplot(data = prds2, aes(x = age, y = Estimate, color = seconds_f)) + 
#   geom_point()

Soil Carbon Example

In the book the use the example of the Lepus and Lynx populations.

The simple differential equation proposed is:

\[ \frac{dH}{dt} = H_t \times (\text{birth rate}) - H_t \times (\text{death rate}) \]

For soil carbon the analogy would be

\[ \frac{dSOC}{dt} = SOC_t \times (\text{assimilation}) - SOC_t \times (\text{decomposition}) \]

SOC-SIMPLE

SOC-CENTURY

SOC-EQ

General Linear Madness

Fernando Miguez